/*
By making sine and cosine table
inline fftType fftSinR(int n, int d){   return fftSin(n*(M_PI_4/d)); }
inline fftType fftCosR(int n, int d){   return fftCos(n*(M_PI_4/d)); }

to evaluate

sin((pi/4 )*(j/w)) or cos((pi/4 )*(j/w)) (w = 2^k, j/w <= 8).

Used in "void rdft(int n, int isgn, fftType *a, int *ip, fftType *w);" only.

Both j and w integers. The maximum value of w = fftMaxArraySize
in "sfft.h)
const uint fftMaxArraySize = 1u << FFT_MAX_SIZE_BITS;
*/
#include "sn.h"

long double ldpow10(const int n);
long double doubleLD(const SDouble& sd, int err);

static fftType  *sine = NULL;   //table of sine
static fftType  *cosine = NULL; //table of cosine

static int L; // FFT_MAX_SIZE_BITS;
static long Nmax; // fftMaxArraySize; Nmax = 1 << L
static long Nmax2;      // Nmax/2
static long Nmax4;      // Nmax/4
static long Nmax3_4;    // 3*Nmax/4
static long Nmax8;      // Nmax/8
static long tableSize = 0;
static void SetNSize(int maxsizebit){
    L = maxsizebit;
    Nmax = 1 << L;  Nmax2 = Nmax/2; Nmax4 = Nmax/4; Nmax3_4 = 3*Nmax4;  Nmax8 = Nmax/8;
    tableSize =Nmax4 +1;
}
long FFTSineTableSize(){
    return tableSize;
}
long FFTSineTableUsedMemory(){
    return tableSize*sizeof(fftType);
}
/*--------------------------------------------------
Here we have the table of sine[k] = sin(2*k*pi/N) where N = 2^L .
We want sin (pi/4)*(j/w) from the table sine[k] (0<= k <= N/4).
Map (pi/4)*(j/w) = 2*k*pi/N, k = ((N/8)*j)/w
(pi/4)*(j/w) = (pi/2)*(j/(2*w))=(pi/2)*n (n = j/(2*w))
----------------------------------------------------*/
static void ErrorEnd(){
    cerr << "out of range in sine table." << endl;
    exit(-1);
}

fftType fftSinR(long j, long w){
    if(!j) return 0.0;
    long k =(Nmax8/w)*j; // (Nmax8*j)/w; pay attention to over flow

    if(j>8*w || k > Nmax || k <0) ErrorEnd();
    if(k <= Nmax4) return sine[k];  // k<N/4
    if(k <= Nmax2) return sine[Nmax2-k];// N/4<k<N/2
    if(k <= Nmax3_4) return -sine[k-Nmax2]; // N/2<k<(3/4)N
    return -sine[Nmax-k];   // (3/4)N<k<N
}

fftType fftCosR(long j, long w){
    if(!j) return 1.0;
    long k = (Nmax8/w)*j; //(Nmax8*j)/w;
    if(j>8*w || k > Nmax || k < 0) ErrorEnd();
    if(k <= Nmax4) return sine[Nmax4-k];// k<N/4
    if(k <= Nmax2) return -sine[k-Nmax4];// N/4<k<N/2
    if(k <= Nmax3_4) return -sine[k-Nmax2];// N/2<k<(3/4)N
    return sine[k-Nmax3_4]; // (3/4)N<k<N
}
fftType fftSin_n(long k){ return sine[k]; } // k < Nmax4
fftType fftCos_n(long k){ return sine[Nmax4-k]; }

/***************************************
It makes the largest sine table first.
Set "maxsizebit" a FFT_MAX_SIZE_BITS.
maxsizebit = 0 : free memory
maxsizebit : same value previous one, do nothing
****************************************/
int MakeSinTable(int maxsizebit){
    if(!maxsizebit){    //free memory
        if(tableSize){  delete[] sine;  sine = NULL;    tableSize = 0;  }
        return -1;
    } else if(L == maxsizebit) return 1;  // same table

    SetNSize(maxsizebit);
    long p;

    delete[] sine;
    delete[] cosine;
    if( (sine = new fftType[tableSize]) == NULL) return 1;
# if 0
    long q,r;
    if( (cosine = new fftType[tableSize]) == NULL ) return 1;
// cr[] and sr[] : table of cos(2*pi/(2^r)) and sin(2*pi/(2^r)) r = 0,...,L
// Almost no this part uses time.
    fftType cr[L+1], sr[L+1];
    cr[0] = 1.0;    cr[1] = -1.0;   cr[2] = 0.0;
    sr[0] = sr[1] = 0.0;            sr[2] = 1.0;
    fftType pi = M_PI, x;
    long m = 4;
    for(r = 3; r < L; r++){
        x= pi/m;
        cr[r] = fftCos(x);
        sr[r] = fftSin(x);
        m *= 2;
    }
    
/***************
 table of sine[n] = sin(2*pi*n/N) and cosine[n] = cos(2*pi*n/N)
  (0 < n < N/4) by an addtion theorem
 Let n = 2^r + q = p + q , p = 2^r
 sine[p+q]
 = sin(2*pi*(p+q)/N)= sin(2*pi*p/N)cos(2*pi*q/N)+cos(2*pi*p/N)sin(2*pi*q/N)
 = sin(2*pi/2^(L-r))*cos(2*pi*q/N) + cos(2*pi/2^(L-r))sin(2*pi*q/N)
 = sr[L-r]*cosine[q]+cr[L-r]*sine[q]
 sine[p]=sin(2*pi*2^r/2^L) = sin(2*pi/2^(L-r)) = sr[L-r]
then
 sine[p+q] = sine[p]*cosine[q]+cos[p]*sine[q]
same way
 cosine[p+q] = cosine[p]*cosine[q]-sine[p]*sine[q];
****************/
Timer timer(true);
    r = 0;  sine[0] = 0.0;      sine[n4] = 1.0;     cosine[0]=1.0;  cosine[n4]=0.0;
    for(p = 1; p < n4 ; p *= 2){
        sine[p] = sr[L -r];     cosine[p] = cr[L -r];       r++;
        for(q = 0; q < p ; q++){
            sine[p+q]   = sine[p]*cosine[q]+cosine[p]*sine[q];
            cosine[p+q] = cosine[p]*cosine[q]-sine[p]*sine[q];
        }
    // here p+q<= N/4   
    }
    delete[] cosine;    cosine = NULL;  // cosine is not necesary
double t = timer.Stop();
cout << "in MakeSinTable() t = " << t << endl;
#else
// sine[n] = sin(2*pi*n/N)
Timer timer(true);
    fftType x, pi2N = M_PI;
    pi2N *= 2; // 2*pi
    pi2N /= Nmax; // 2*pi/N
    sine[0] = 0.0;      sine[Nmax4] = 1.0;
    for(p = 1; p < Nmax4 ; p++){
        x = p*pi2N;
        sine[p] = fftSin(x);
    }
double t = timer.Stop();
cout << "in MakeSinTable() t = " << t << endl;
    return 0; // OK
}
#endif

/*----------------------- end of sine table---------------------------*/